*****************************************************************
*** Do file to replicate the analyses for the three illustrations in Laron K Williams, "Compression, Temporal Dependence, and the Sensitivity of Quantities of Interest"
***
*** Created: 7-29-14
*** Modified: 3-6-18
***
*****************************************************************

*** Set up the working directory
*cd ""

*****************************************************************
*** Ways and Weeks (2014, AJPS)
*****************************************************************
use "Way and Weeks 2014\WayWeeksAJPS.dta", clear

*** Table 1 (page 714)
btscs pursueonly year ccode, g(time)
gen time2=time*time
gen time3=time2*time

*** Linear probability model
xtreg pursueonly persdumjlw_lag land gdppercap time time2 time3, re
di _b[persdumjlw_lag]	

*** Logit model
xtlogit pursueonly persdumjlw_lag land gdppercap time time2 time3, nolog
keep if e(sample)

*** Predict
predict base, pr

tempvar pr_0 pr_1 
replace persdumjlw_lag = 0
predict `pr_0', pu0

replace persdumjlw_lag = 1
predict `pr_1', pu0
gen pe = `pr_1' - `pr_0'

qui sum pe
local ape = round(r(mean), 0.001)
local sdpe = round(r(sd), 0.001)

di "Average partial effect = " `ape'
di "Standard deviation of partial effect = " `sdpe'

sum base `pr_0' `pr_1' pe
	
*** QIs across t	
bys time: egen ape_t = mean(pe)
bys time: egen sdpe_t = sd(pe)
bys time: egen minpe_t = min(pe)
bys time: egen maxpe_t = max(pe)

twoway (scatter pe time, jitter(3)) (line ape_t time), ytitle("{&Delta}Pr(Y=1)") xtitle("Years Since Last Pursuit") legend(off) 
twoway (scatter pe time, jitter(3)) (line sdpe_t time)

*** Create a dataset to generate the figure in R
preserve
	keep pe time ape_t sdpe_t minpe_t maxpe_t
	saveold "Way and Weeks 2014\Data\WW2.dta", replace version(12)
	
	duplicates drop time, force
	list
restore


*****************************************************************
*** Cunningham (2013, AJPS), Model 2
*****************************************************************
use "Cunningham 2013\CunninghamAJPS.dta", clear
gen id = _n

*** Logit replication
logit acdcivilwar1 logfactions prevconcessions_l democracy kin yrsnocivilwar _spline1_cw _spline2_cw _spline3_cw , robust cluster(kgcid)
keep if e(sample)

*** Linear probability model
reg acdcivilwar1 logfactions prevconcessions_l democracy kin yrsnocivilwar _spline1_cw _spline2_cw _spline3_cw , robust cluster(kgcid)
di _b[logfactions] * 3.66

*** Replication of substantive effects
estsimp logit acdcivilwar1 logfactions prevconcessions_l democracy kin yrsnocivilwar _spline1_cw _spline2_cw _spline3_cw , robust cluster(kgcid)

* Logged Factions
setx mean
setx prevconcessions_l 0 democracy 0 kin 0 
simqi, fd(prval(1)) changex(logfactions 0 3.66)

* Democracy
setx mean
setx prevconcessions_l 0 democracy 0 kin 0 
simqi, fd(prval(1)) changex(democracy 0 1)

*** Overall average partial effect and standard deviation of partial effect
logit acdcivilwar1 logfactions prevconcessions_l democracy kin yrsnocivilwar _spline1_cw _spline2_cw _spline3_cw , robust cluster(kgcid)

predict base, pr

tempvar pr_0 pr_1 
replace logfactions = 0
predict `pr_0', pr

replace logfactions = 3.66
predict `pr_1', pr
gen pe = `pr_1' - `pr_0'

qui sum pe
local ape = round(r(mean), 0.001)
local sdpe = round(r(sd), 0.001)

di "Average partial effect = " `ape'
di "Standard deviation of partial effect = " `sdpe'

sum base `pr_0' `pr_1' pe

preserve
	keep id pe
	rename pe pe_lf
	sort id
	tempfile data
	save `data', replace
restore	
	
bys yrsnocivilwar: egen ape_t = mean(pe)
bys yrsnocivilwar: egen sdpe_t = sd(pe)
bys yrsnocivilwar: egen minpe_t = min(pe)
bys yrsnocivilwar: egen maxpe_t = max(pe)

twoway (scatter pe yrsnocivilwar, jitter(3)) (line ape_t yrsnocivilwar), ytitle("{&Delta}Pr(Y=1)") xtitle("Time Since Civil War Incidence") legend(off) 
twoway (scatter pe yrsnocivilwar, jitter(3)) (line sdpe_t yrsnocivilwar)

preserve
	rename yrsnocivilwar time
	keep pe time ape_t sdpe_t minpe_t maxpe_t
	saveold "Cunningham 2013\Data\Cunningham.dta", replace version(12)
	
	duplicates drop time, force
	list time *_t
restore

preserve
	rename yrsnocivilwar time
	keep time 
	saveold "Cunningham 2013\Data\chist.dta", replace version(12)
restore

*** Average partial effects across other variables
foreach v of varlist prevconcessions_l democracy kin {
	bys `v': egen ape_`v' = mean(pe)
	bys `v': egen sdpe_`v' = sd(pe)
	bys `v': egen minpe_`v' = min(pe)
	bys `v': egen maxpe_`v' = max(pe)
	
	preserve
		keep `v' *_`v'
		duplicates drop `v', force
		list `v' *_`v'
	restore
}


*** Average partial effect of democracy
use "Cunningham 2013\CunninghamAJPS.dta", clear
gen id = _n

preserve
	qui logit acdcivilwar1 logfactions prevconcessions_l democracy kin yrsnocivilwar _spline1_cw _spline2_cw _spline3_cw , robust cluster(kgcid)
	keep if e(sample)

	tempvar pr_0 pr_1 
	replace democracy = 0
	predict `pr_0', pr

	replace democracy = 1
	predict `pr_1', pr
	gen pe = `pr_1' - `pr_0'

	qui sum pe
	local ape = round(r(mean), 0.001)
	local sdpe = round(r(sd), 0.001)

	di "Average partial effect = " `ape'
	di "Standard deviation of partial effect = " `sdpe'

	sum `pr_0' `pr_1' pe
restore

*****************************************************************
*** Flores-Macias and Kreps (2013, APSR)
*****************************************************************
use "Flores-Macias and Kreps 2013\Flores-MaciasKrepsAPSR.dta", clear

*** Linear probability model
preserve
	reg taxdummy party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3 if severity>3, vce(robust)
	di _b[party]
restore

*** Flores-Macias and Kreps replication
preserve
	* Table 3 (page 842)
	estsimp logit taxdummy party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3 if severity>3, vce(robust)

	* Set all continuous variables to means and dummy variables to medians 
	setx divided median
	setx lagdebtgdp mean
	setx laggdpgrowth mean
	setx lagPerChangeDefExp mean
	setx type median
	setx yearsnotax mean
	setx cubicspline1 mean
	setx cubicspline2 mean
	setx cubicspline3 mean
	setx electionyear median
	setx party median
	
	simqi, listx
	
	* Calculate the marginal effect of the following variables on the probability of a war tax 
	simqi, fd(pr) changex(party 0 1)
	simqi, fd(pr) changex(severity p25 p75)
	simqi, fd(pr) changex(mid_sidea 0 1)
	simqi, fd(pr) changex(lagPerChangeDefExp p25 p75)
	simqi, fd(pr) changex(electionyear 0 1)
	simqi, fd(pr) changex(laggdpgrowth p25 p75)
	simqi, fd(pr) changex(type 0 1)
	simqi, fd(pr) changex(divided 0 1)
	simqi, fd(pr) changex(laginflationrate p25 p75)
	simqi, fd(pr) changex(lagdebtgdp p25 p75)
restore

qui logit taxdummy party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3 if severity>3, vce(robust)
keep if e(sample)
local end = e(N)
mat bl = e(b)'

*sum party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3, det

preserve
	duplicates drop yearsnotax, force
	keep yearsnotax cubic*
	sort yearsnotax
	mkmat yearsnotax cubicspline1-cubicspline3, mat(S)
	local r = rowsof(S)
	local rminusone = `r' - 1
restore

mat list S

di "Average case approach has the following values for the simulation scenario:"
di "Years no tax = 14.4"
di "Spline #1 = 5.8"
di "Spline #2 = 3.7"
di "Spline #3 = 1.4"

estsimp logit taxdummy party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3 if severity>3, vce(robust)

* Not quite right
setx mean
setx divided 0 electionyear 0 mid_sidea 1 type 0 severity 4 
simqi, fd(prval(1)) changex(party 0 1)

qui foreach i in fd fd_lo fd_hi xaxis {
	gen `i' = .
}

* Loop across (correct) values of t
qui foreach s of numlist 1(1)`r' {
	tempvar fd_`s'
	local s1 = S[`s',1]
	local s2 = S[`s',2]
	local s3 = S[`s',3]
	local s4 = S[`s',4]

	setx mean
	setx divided 0 electionyear 0 mid_sidea 1 type 0 severity 4  yearsnotax `s1' cubicspline1 `s2' cubicspline2 `s3' cubicspline3 `s4'
	simqi, fd(genpr(`fd_`s'') prval(1)) changex(party 0 1)
	sum `fd_`s'', meanonly
	replace fd = r(mean) in `s'
	
	_pctile `fd_`s'', p(2.5 97.5)
	replace fd_lo = r(r1) in `s'
	replace fd_hi = r(r2) in `s'
	
	replace xaxis = `s1' in `s'
}

preserve
	keep yearsnotax
	saveold "Flores-Macias and Kreps 2013\Data\FMKdist.dta", replace version(12)
restore

*** Observed-values approach
preserve
	qui logit taxdummy party divided electionyear mid_sidea lagdebtgdp laginflationrate laggdpgrowth lagPerChangeDefExp type severity yearsnotax cubicspline1 cubicspline2 cubicspline3 if severity>3, vce(robust)
	keep if e(sample)

	predict base, pr
	
	tempvar pr_0 pr_1 
	replace party = 0
	predict `pr_0', pr

	replace party = 1
	predict `pr_1', pr
	gen pe = `pr_1' - `pr_0'

	qui sum pe
	local ape = round(r(mean), 0.001)
	local sdpe = round(r(sd), 0.001)

	di "Average partial effect = " `ape'
	di "Standard deviation of partial effect = " `sdpe'

	sum base `pr_0' `pr_1' pe
	sum pe, det
	hist pe
	
	gen pe_fmk = .061
	bys yearsnotax: egen pe_mn = mean(pe)
	keep pe* yearsnotax 
	saveold "Flores-Macias and Kreps 2013\Data\Observed-Case Partial Effects--FMK.dta", replace version(12)	
restore

*** They underestimate the true effect by using the average-case approach!
